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ON RECONSTRUCTING TRAJECTORIES IN 
THE VENUS LOWER ATMOSPHERE 


P. Argentiero 
G. Wyatt 


ABSTRACT 

This paper utilizes a Monte Carlo technique in 
order to demonstrate the feasibility of processing in- 
situ measurements of temperature, pressure, and 
molecular weight in order to reconstruct trajectories 
in the Venus lower atmosphere. The technique as- 
sumes that the Venus lower atmosphere obeys the 
ideal gas law and the hydrostatic equation. Time cor- 
relations in the data are assumed to exist. It is also 
shown that the errors in trajectory reconstruction 
are due mostly to noise on the data rather than inac- 
curacies in the numerical technique. 
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ON RECONSTRUCTING TRAJECTORIES IN 
THE VENUS LOWER ATMOSPHERE 


INTRODUCTION 

The problem of reconstructing a trajectory in an atmosphere has been dis- 
cussed in the literature. See for instance Preslin [1], Peterson [2], Sommer 
and Boissevain [3], and Seiff [4]. These studies assume that in situ measure- 
ments of accelerations temperature and pressure are available. In the Venus 
lower atmosphere the density is so great that measurements from on board ac- 
celerations are not likely to be useful. Thus trajectory reconstruction in the 
Venus lower atmosphere must rely on range rate measurements and in-situ 
measurements of pressure, temperature, and molecular weight. The use of 
these data types in a standard parameter estimation mode requires that the 
planetary atmosphere be characterized by a relatively small number of 
parameters. This is usually done by assuming that the atmosphere satisfies 
the hydrostatic equation and the ideal gas law with the temperature profile piece- 
wise linear with a small number of breakpoints. In this case the surface pres- 
sure and the temperature at the breakpoints along with the height of the break- 
points are parameters which define an atmosphere. A standard Kalman filter 
approach can then be implemented and an estimate of the trajectory can be ob- 
tained. The atmospheric parameters can be estimated as part of the state or 
they can be placed in a consider mode where they are not estimated but their 
associated uncertainties are permitted to have an influence on the filter behavior. 
The ensemble properties of such a filter can be studied in an error analysis 
mode and if the modeling assumptions are valid and if the non-linearities of the 
problem are not severe, the promulgated covariance matrix of such an error 
analysis represents an accurate measure of the statistical quality of the estimate. 

Unfortunately there is no reason to believe that the non linearities of the 
problem are negligible. There are also present some obvious modeling errors 
in this approach. For instance it is doubtful if the temperature profile of the 
lower Venus atmosphere is approximately piecewise linear. The significance of 
such modeling errors is not now known. Consequently the results of an ensemble 
calculation associated with an error analysis of the parameter estimation tech- 
nique is not to be trusted as an accurate indicator of the statistical quality of the 
technique. Another method for obtaining information concerning the statistical 
spread of the parameter estimation procedure is to perform a Monte Carlo study. 
But the computational algorithms for the implementation of the parameter esti- 
mation technique are too lengthy to permit such a possibility. One is forced to 
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conclude then that there is no effective way of determining the statistical quality 
of the standard parameter estimation technique as applied to trajectory determin- 
ation in the Venus lower atmosphere. 

There is a way out of this impasse. If one is willing to forego the direct 
use of range rate data, then a computational algorithm can be devised which 
utilizes temperature, pressure, and molecular weight data in order to estimate 
the trajectory. This algorithm is sufficiently simple that its ensemble proper- 
ties can be accurately determined by means of a Monte Carlo study. This algo- 
rithm is obtained by assuming that measurements are unbiased and that the 
atmosphere obeys the hydrostatic equation and the ideal gas law. No assump- 
tions concerning the temperature profile are necessary. This paper is a report 
on a Monte Carlo study of such an estimation procedure. 


ESTIMATING ATMOSPHERIC TRAJECTORIES WITH PRESSURE, 
TEMPERATURE, AND MOLECULAR WEIGHT DATA 

Assume that an atmosphere is perfectly mixed, and that it obeys the per- 
fect gas law and the hydrostatic equation. Let P(H), T(H), p (H) be respectively 
the pressure, temperature and density of the atmosphere at a height H above 
the surface. Then 


and 


p (H) R T(H) 

P(H) = „ 

M 


dP(H) _ 
dH 




( 1 ) 

( 2 ) 


where R is the universal gas constant, M is the molecular mass of the gas and g 
is the gravitational constant of the planet. (For present purposes the variation 
of g with H can be ignored.) 

Assuming a ninety degree flight path angle, the trajectory of a probe in the 
Venus lower atmosphere can be described by a height versus time function H(t). 
Define the following functions 


P(t) = P CH( t )] 

(a) 


T(t) = T [H(t)] 

(b) 

(3) 

p(t) =yofe(t)] 

( c ) 
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Equations 1 and 2 yield 


dH(t) _ -RT(t) 
dP(t) " gP(t)M 


By multiplying both sides of equation 4 by DP(t)/Dt one obtains 


and equation 5 yields 


dH(t) _ -R T(t) dP(t ) 
dt " gM P(t) dt 


H(t) 


JL r to T(t) dP(t) dt 

gM J P(t) dt 


(5) 


(6) 


where t 0 is the time of impact. Equation 6 provides an effective algorithm for 
reconstructing the height of the probe as a function of time. 

The salient difficulty with the approach suggested by equation 6 is that the 
readings of temperature and pressure are available in discrete rather than 
continuous form. Thus both the derivative of pressure with respect to time and 
the integration of equation 6 must be obtained according to some numerical 
procedure. It is by no means obvious which numerical techniques would yield 
the smallest errors in the reconstruction of the trajectory. The difficulty is 
further compounded by the fact that there is noise on the data. Thus the error 
in trajectory reconstruction is ultimately the result of the interaction between 
the deterministic features of the numerical technique and the ensemble features 
of the data noise. The only systematic way to decide which of a set of numerical 
techniques yields the smallest errors in trajectory reconstruction is to test each 
in a Monte Carlo mode and to choose the numerical technique with the smallest 
critical error values as the most accurate. In the process one not only de- 
termines the best numerical technique but he also establishes the critical values 
of the error distribution of the resultant trajectory determination if the technique 
were utilized. 


THE MONTE CARLO SIMULATION 

In order to perform a Monte Carlo simulation, one must first construct a 
deterministic model of reality. In this case reality consists of a model of the 
Venus atmosphere and a probe with a certain trajectory which is sampling the 
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temperature, pressure, and molecular weight of the atmosphere according to a 
time schedule. 

Since the atmosphere is assumed known and the trajectory and data sampling 
schedule are fixed, one knows the true values of the collected data. A random 
number generator is used to add noise to the true values of the data. This noise 
is so constructed as to follow a predetermined probability law. After the noise 
has been added, the resultant numbers are processed by the computational 
algorithm in use and an estimate of the trajectory is obtained. Since the true 
value of the trajectory is assumed known, one can' determine the errors which 
the algorithm yields. The entire process can be repeated many times and an 
estimate of, say, the 95% critical value of the error in trajectory estimation at 
any given height can be obtained. The noise on the measurements was assumed 
to be zero mean and normally distributed. The assumed standard deviations as 
obtained from John Ainsworth of the Laboratory For Planetary Atmospheres, 
G.S.F.C., are displayed in Tables I and n. 


Table I 

Standard Deviation of Pressure Measurements 


Range in Earth 

Standard Deviation as Percent of 

Atmosphere 

Largest Value in Interval 

90 - 20 

.1% 

20 - 4 

.2% 

4 - .08 

A% 

.08 - .05 

.5% 


Table n 

Standard Deviation of Temperature Measurements 


Range in Kelvin 

Standard Deviation as Percent of 
Length of Interval 

800 - 600 K 

1% 

600 - 400 K 

.75% 

400 - 300 K 

.60% 

300 - 200 K 

.5% 


The G.S.F.C. 3609 model for the Venus atmosphere was used for our simu- 
lations. The probe trajectory was that of an entry probe with a 90° flight path 
angle and a parachute opening at 50 km above the surface. The reconstruction 
process is assumed to begin at 75 km above the surface. The trajectory was 
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constructed in the G.S.F.C. 3609 atmosphere by using the planetary explorer 
main probe configuration and assuming the quasi-static approximation. This 
approximation implies that the downward force of gravity is always in approxi- 
mate equilibrium with the upward drag force. In effect this provides the velocity 
of the probe as a function of the density of atmosphere and makes for an easy 
trajectory reconstruction. Recently, this technique has been compared to a 
fourth order Runge Kutta technique and has been found to be valid in the Venus 
lower atmosphere, Reference [5]. 

Pressure, temperature, and molecular weight were assumed to be sampled 
a total of seventy-six times. The time intervals between measurements were 
constant and the first measurements were taken at 75 km and the last at the 
surface. Certain systematic errors in the measurement processes induce time 
correlations in the data. To model these measurements as independent is to 
make an optimistic error since this exaggerates the information content of the 
data. In this study the authors assume a positive correlation of .8 between a 
temperature or pressure measurement and the previous temperature or pres- 
sure measurement. These correlations become respectively .6, .4 and .2 when 
the time lag becomes 2, 3 and 4. When the time lag is 5 or greater we assume 
the measurements are independent. No cross correlation between the tempera- 
ture and the pressure measurements are assumed to exist. The standard devia- 
tion of the noise on the measurements of molecular mass is assumed to be .5% 
of the nominal value of 43 as used in the 3609 model of the Venus atmosphere. 

A systematic error of .001 cm/sec 2 was assumed to have been made in the cal- 
culation of the surface gravitational acceleration of Venus. 

The standard way to obtain Monte Carlo estimates of 95% critical values of 
errors at various heights is to execute the algorithm under investigation N times 
with the noise on the data chosen each time according to the probability law 
defined above. Each of the N trajectories so obtained can be compared with the 
true trajectory and at any given height, the 95% largest error from the ensemble 
of N errors is the optimal estimate of the 95% critical value of the error of the 
trajectory at that height. Unfortunately this optimal Monte Carlo technique when 
applied to this problem, leads to intractible problems both with regard to com- 
puter time and core storage. Thus a sub optimal Monte Carlo technique was 
utilized. The details of this technique are discussed in the appendix. This sub- 
optimal technique does not sacrifice rigor and one can be sure to a confidence 
level of 99% that the estimate of a critical value it yields is larger than the 
correct critical value. In this sense the technique can be used to produce con- 
servative estimates of critical values. 
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RESULTS AND CONCLUSIONS 

Figures 1 through 4 provide graphs of respectively the 90%, 95%, 97% and 
99% critical values of the error in trajectory reconstruction as a function of 
height above the surface. The algorithm used was implied by equation 6 with 
the derivative of pressure with respect to time calculated by assuming that 
pressure is exponential with time and with the integration implemented with 
the aid of a Simpson's rule algorithm. When the pressure derivatives are cal- 
culated by means of a simpler linear approximation, the errors increase by al- 
most an order of magnitude. The use of a Simpson's rule in order to implement 
the integration instead of a simpler trapezoidal rule purchased an improvement 
in the errors of about 5%. Obviously the way the pressure derivatives are esti- 
mated is far more significant in the trajectory reconstruction than how the 
integration is performed „ 

From Figures 1 through 4 it is evident that at least to a height of 50 km 
above surface the processing of temperature and pressure data according to 
equation 6 yields a good estimate of the probe's trajectory. Above this height, 
the modeling assumptions, primarily the assumption that pressure measure- 
ments are exponential in time, become increasingly invalid. This above- 
mentioned exponential assumption rests on the assumptions that pressure is 
exponential with height and that probe velocity is constant. Clearly the greater 
the height above surface the greater is the acceleration the probe processes and 
the more questionable the assumption of constant velocity becomes. This is the 
primary impediment to the extension of this trajectory estimation technique to 
greater heights above surface. Of course, the technique's effectiveness could 
be extended considerably if the data acquisition rates were altered so that the 
higher the probe is above surface the greater the data acquisition rate. Although 
the authors have not done so, the possibility of extending this technique to much 
greater heights by means of a varied data acquisition rate could be systematically 
studied with the aid of the Monte Carlo techniques demonstrated in this report. 

It is of some value to understand how much of the errors in trajectory re- 
construction is due to noise on the data and how much is due to inaccuracies 
inherent in the numerical technique that one decides to use. In the present case 
such information is easy to obtain. Figure 5 provides a graph of the error ob- 
tained in trajectory reconstruction with perfect data. In this case the errors are 
due entirely to inaccuracies in the numerical technique. It is evident in compar- 
ing the information in Figure 5 with the information in any of the other figures 
that noise on the date rather than errors in the numerical technique is primarily 
responsible for errors in trajectory reconstruction. For instance compare the 
95% critical error value in Figure 2 at 40 km, to the deterministic error at 40 
km as given in Figure 5. The value from Figure 5 is smaller than the corre- 
sponding value from Figure 2 by an order of magnitude. The values of Figure 5 
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Figure 2. 95% Critical Error Value in Trajectory Reconstruction 
as a Function of Height Above Surface 
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Figure 3. 97% Critical Error Value in Trajectory Reconstruction 
as a Function of Height Above Surface 
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Figure 5. Deterministic Error Value in Trajectory Reconstruction 
as a Function of Height Above Surface 
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ma y be thought of as the values toward which the errors in trajectory recon- 
struction would asymptotically tend as the data became increasingly more ac- 
curate. As such these values represent the limit of how much can be purchased 
by improvement in data quality. Clearly there is substantial room for 
improvement. 
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APPENDIX 


OPTIMAL AND SUB OPTIMAL MONTE CARLO TECHNIQUES 
FOR ESTIMATING CRITICAL VALUES 


One usually employs a Monte Carlo technique in order to estimate a 
parameter associated with a probability distribution. In this study we have 
been concerned with the Monte Carlo estimation of critical values of a distri- 
bution. What follows below is an analysis of both an optimal and a sub optimal 
Monte Carlo technique for estimating critical values. 

Given a sample of N independent valuations of a random variable, the 
"optimal" processing of this data in order to estimate the critical value of the 
random variable which corresponds to probability P can be very simply 
described. Let K be the smallest integer which is larger than NP. Linearly 
order the samples in descending order of magnitude. The K number in such 
a list is taken as the estimate of the critical value corresponding to probability 
P. We have yet to define in what sense such an estimate is optimal nor have we 
specified how the quality of the estimate depends on the sample size N. 

In order to discuss these questions, it is first necessary to analyze the proba- 
bility structure of the following experiment. Choose a probability P according 
to a uniform probability law. Then perform N independent Bernoulli trials with 
probability of success P. Divide the number of successes in N trials by N and 
call the result P'. The result of the experiment is the two tuple (P, P'). If N 
is sufficiently large, then an application of the normal approximation to the bi- 
nomial distribution yields the conditional probability density function (p.d.f.) of 
P given that P* was fixed at some value a , 0 < a < 1 as 


F(P) 



exp 


1 

2 


P - a “I 



2 


(1) 


Whenever one performs N Bernoulli trials where the probability of success is 
unknown, and the ratio of the number of successes to N is a, then equation 1 may 
be considered as an approximation to the p.d.f. of the true probability of 
success P. 
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With regard to the previously defined estimation procedure, the Bernoulli 
trials may be thought of as the independent sampling from the distribution and 
success in a trial can be defined as the probability of exceeding our estimated 
critical value. Since we choose this critical value such that the ratio of suc- 
cesses to trials is approximately a, equation 1 is the approximate probability 
density function of the true probability of the random variable in question ex- 
ceeding our estimated critical value. The maximum value of the right side of 
equation 1 is achieved when P = a. Thus the maximum likelihood estimate of 
the probability level associated with the estimated critical value is a. It is in 
this sense that the estimate is optimal. 

We would like to obtain a region about the point P = a which would contain 
say 95% of the area under the curve defined by equation 1. If such regions could 
be constructed, then questions concerning the accuracy of the estimation tech- 
nique could be adequately answered. To do this we substitute a probability 
density function whose integrated values are tabulated and which has the prop- 
erty that a 95% critical region about a contains more than 95% of the area under 
the p.d.f. defined byequation 1. A hint as to how this could be done is obtained 
by observing that the following inequality is valid for 0 < P < 1 

(2) 

2 /N f N 


The substitution of the left side of inequality 2 for the expression on the right 
side transforms the right side of equation 1 into 


g(P) 


1 


v/277 —L= 

2 /N 


exp 


1 

2 



(3) 


Equation 3 provides the p.d.f. of a normal random variable with expectation 
a and standard deviation equal to l/2 /N. In order to determine if this distri- 
bution has the necessary properties as a bound on the p.d.f. of equation 1, we 
must discover for what range of values g(P) bounds F(P). This is done by ob- 
taining an expression in P for the log of the ratio of g(P) to F(P). The values of 
P for which this expression is positive determines the region for which g(P) 
bounds F(P). After considerable algebraic manipulation, the inequality to be 
solved takes the following form. 


0 < P(1 - P) LN 2 /P(l - P) + (P - a ) 2 N (4P 2 - 4P + 1) (4) 


14 



It is readily shown that if P is in the region 



( 5 ) 


then inequality 4 is satisfied. The area under g(P) in this region is .05. The 
area under F(P) in this region must be less than this value. Thus the p.d.f. of 
equation 3 has the desired properties as a bound. 

It now becomes easy to determine the relationship between the accuracy of 
the estimation technique and the sample size N. If a is the estimated probability 
level associated with the critical value s then one is certain to a confidence level 
of .95 that the true probability level will not differ from a by more than l// N. 

To be certain to this confidence level that the estimated probability level is not 
incorrect by more than .01, one needs a sample size of 10,000. If a confidence 
level of .99 is desired, then a sample size of 15,000 is needed. 

The previous discussion has outlined the standard Monte Carlo technique 
for obtaining optimal estimates and has indicated the relationship between 
sample size and the accuracy. The technique has shortcomings, one practical 
and the other theoretical. To start with the practical difficulty, the technique 
can be demanding computationally if large sample sizes are involved. All the 
sample values must be stored in the computer at one time in order to implement 
any ofthe algorithms whose execution times are not prohibitive and in some 
circumstances this is not possible. The theoretical difficulty is that the optimal 
technique permits us to make probabilistic statements concerning the true 
probability level associated with our estimate of, say, the critical value associ- 
ated with the .05 probability level. But we can say nothing probabilistically about 
the correct value of the critical value associated with the .05 probability level. 
For instance it could be useful in many circumstances to specify an interval 
and be certain to a given confidence level that the true critical value is in the 
interval. There exist so-called sub optimal Monte Carlo estimation techniques 
which can provide such information and which algorithmically are more tract- 
ible than the optimal technique. One such technique, the one that was utilized 
to obtain the results quoted in the body of this report will now be described. 

Suppose one obtains K batches of samples each containing N values and 
all the values are assumed to be independent. For each of the K batches we can 
process N samples in the batch in the optimal fashion described previously and 
thus obtain K numbers X ; , i ^ K all of which are estimates of the critical value 
X associated with a probability level a. Each of these estimates is equally 
likely to be less than as greater than the correct critical value X. Thus the 
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number of our K estimates which exceed X is a binomial random variable with 
expectation K/2 and variance K/4. Using the usual normal approximation to the 
binomial distribution we can infer to a confidence level of .99, that X is smaller 
than at least M 2 of the K values where M 2 is the largest integer which is 
smaller than K/2 + 3 / K/2 and that X is larger than of the values where M x 
is the smallest integer which is larger than K/2 - 3 v K/2. The significance of 
these facts becomes apparent when we state them in the following logically 
equivalent form. Suppose we order the K estimates of X in ascending order of 
magnitude. Then one is sure to a confidence level of .99 that the correct criti- 
cal value X lies between the value in the sequence indexed by M ( and the value 
in the sequence indexed by M 2 . Thus by this technique we have succeeded in 
structuring an interval in which we are certain to a confidence level of .99 the 
correct critical value X resides. The algorithmic advantages of this technique 
rest chiefly in the fact that only one batch of values at a time need be stored in 
the computer. There are also advantages in terms of computing speed. 

While rigorous statistical results concerning the accuracy of this technique 
are difficult to obtain, it is intuitively evident how the batch size N and the num- 
ber of batches K influence the accuracy. The number N determines the accuracy 
of each individual estimate, that is, N determines how widely scattered the esti- 
mates are about the true value X and thus how large the interval which is ex- 
pected to contain X will be. If N is small these estimates will be widely scattered 
and the interval will be much larger than necessary and thus too conservative. 

The impact of the number of batches K on the accuracy of the technique is more 
subtle but nevertheless quite significant. First, if K is too small the normal ap- 
proximation to the binomial distribution becomes invalid and serious errors can 
result. In addition to this factor there is the fact that we are not at liberty to 
choose any confidence level we wish since there are only a finite number of 
intervals that this technique permits one to choose. We pick the interval whose 
confidence level bounds most closely the desired confidence level. The larger 
the value of K which is chosen, the better the quality of this bound. 
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